function IRF_BLP_all = IRF_BLP_estimation(data,lag_type,IRF_hor,Y_var_lagsnon1,Y_var_lags,Y_var_ad,hyperPars,hyperPriorsOptions,scale)
%Created by Viet Hoang Dinh

%----------------- Bayesian LP------------------
%----- For h=0 and 1-----
y=data;
[T,n]=size(y);
nL=lag_type;
nP=lag_type;
nH=IRF_hor;

%-VAR PRIORS-
Ylag=Y_var_lagsnon1;
YprojSet = Y_var_lags; 
Y=Y_var_ad;
nT=size(Ylag,1); 
Y_init=mean(y(1:nL,:)); %average of initial observations

%set prior residual variance (Sigma) using univariate AR(1) residuals
sigmaj=NaN(n,1); 
for j=1:n    
    sigmaj(j)=std( Y(:,j)-[ones(nT,1) Ylag(:,j)]*...
        ([ones(nT,1) Ylag(:,j)]\Y(:,j)) );    
end

%Prior for VAR coefficients: RW
B_init_RW         =zeros(n*nL+1,n);
isRandomWalk      =ones(1,n);
B_init_RW(2:n+1,:)=diag(1.*isRandomWalk); 

%Prior for VAR coefficients: OLS
B_init_OLS=(YprojSet'*YprojSet)\(YprojSet'*Y);

% * * * * * * * * * * * get optimal hyperparameters * * * * * * * * * * * %
%parsAtMode=maxMLikelihoodVAR(Y,YprojSet,B_init_RW,sigmaj.^2,Y_init,hyperPriorsOptions); %Use RW Prior
parsAtMode=maxMLikelihoodVAR(Y,YprojSet,B_init_OLS,sigmaj.^2,Y_init,hyperPriorsOptions); %Use OLS Prior
lambda=parsAtMode.postmax.lambda; % Optimal lambda
lambdaC=hyperPars.lambdaC;

%shrinkage over horizons
optimalLambda=NaN(nH,1); 
optimalLambda(1)=lambda;

%prior variance (coefficients)
Omega_init=inv(blkdiag(1/lambdaC,kron(diag(1:nL).^2,diag(sigmaj.^2))/lambda^2)); 

% Weighting matrices: Bayesian LP
M_matrix=YprojSet'*YprojSet*(Omega_init);
n_coeff=size(M_matrix,1);
Q_matrix=inv(eye(n_coeff)+inv(M_matrix));
Q_matrix_c=eye(n_coeff)-Q_matrix;

% LP and VAR
B_LP=(YprojSet'*YprojSet)\(YprojSet'*Y);
B_VAR_RW=B_init_RW;
B_VAR_OLS=B_init_OLS;

% BLP coefficients matrix
%B_endBayesian=Q_matrix'*B_LP+Q_matrix_c'*B_VAR_RW;  %Use RW Prior
B_endBayesian=Q_matrix'*B_LP+Q_matrix_c'*B_VAR_OLS;  %Use OLS Prior


% Contemporaneous matrix 
if scale==0
    Bzero_Bayesian=chol(cov(Y-YprojSet*B_endBayesian),"lower"); 
elseif scale==1
    Bzero_Bayesian=chol(cov(Y-YprojSet*B_endBayesian),"lower"); 
    for i_Bzero=1:n
    Bzero_Bayesian(:,i_Bzero)=Bzero_Bayesian(:,i_Bzero)./Bzero_Bayesian(i_Bzero,i_Bzero);
    end
end

%Interested coefficients
B_endBLP_int=B_endBayesian(2:end,:);

% Companion matrix
F_varBLP = [B_endBLP_int';eye(n*(nL-1)) zeros(n*(nL-1),n)];
bigeye= [eye(n) zeros(n,n*(nL-1))];

% IRFs at h=0,1;
IRF_BLP_all=zeros(n^2,IRF_hor);


for j_hor= 1:2
    if j_hor==1
    IRF_BLP_all(:,j_hor)= reshape(Bzero_Bayesian,[],1); 
  
    elseif j_hor==2
    IRF_BLP_all(:,j_hor) = reshape(bigeye*F_varBLP*bigeye'*Bzero_Bayesian,[],1);
    
    end    
 
end


%----- For h>1-----
BVARcoeffs_Bayesian=B_endBayesian;
x=[zeros(nL, size(Y,2)); Y];


for h=2:nH-1
    XhprojSet=YprojSet(1:end-h+1,:);
    Xh=Y(h:end,:);
    nT=size(Xh,1);       
    Xh_init=mean(x(1:nP,:)); 
    
 %Use univariate local projection to initialize scale (NW corrected)
    gammaU=NaN(n,1);    
    for k=1:n        
        projCoeffs=XhprojSet(:,[1 k+1:n:(n*nP+1)])\Xh(:,k);

        %univariate projection residuals
        u=Xh(:,k)-XhprojSet(:,[1 k+1:n:(n*nP+1)])*projCoeffs;
        u=bsxfun(@minus,u,mean(u)); 
        GammaU=(u'*u)/nT; %HAC error (T*)covariance estimator        
        nwLags=h-1; %u_{t+h|t} is MA(h-1)
        
        %HAC correction: prior scale
        nwWeights=(nwLags+1-(1:nwLags))./(nwLags+1);
        for j=1:nwLags
            gammaj=(u(j+1:nT,:)'*u(1:nT-j,:))/(nT-j);            
            GammaU=GammaU+nwWeights(j)*(gammaj+gammaj');
        end
        gammaU(k)=sqrt(GammaU); %scalar        
    end
   
%Centered on relevant power of VAR coefficients
  Bh_init_Bayesian=setPriorMean_VAR2(BVARcoeffs_Bayesian(2:end,:),h,n,nL,nP);
  nBh=numel(Bh_init_Bayesian);

 % * * * * * * * * * * get optimal hyperparameters * * * * * * * * * * %
    parsAtMode=maxMLikelihoodBLP(Xh,XhprojSet,Bh_init_Bayesian,gammaU.^2,Xh_init,hyperPriorsOptions,nP,h);
    lambdaP=parsAtMode.postmax.lambda;  %Optimal lambda
    optimalLambda(h)=lambdaP;

 %prior variance
    OmegaH_init=inv(blkdiag(1/lambdaC,kron(eye(nP),diag(gammaU.^2))/(lambdaP^2))); %think of it as inv(Xd'Xd); Xd dummy observations

% Weighting matrices: Bayesian LP
Mh_matrix=XhprojSet'*XhprojSet*(OmegaH_init);
Qh_matrix=inv(eye(n_coeff)+inv(Mh_matrix));
[msg, id] = lastwarn;
warning('off', id);
Qh_matrix_c=eye(n_coeff)-Qh_matrix;

% LP and VAR
Bh_LP=(XhprojSet'*XhprojSet)\(XhprojSet'*Xh);
Bh_VAR_Bayesian=Bh_init_Bayesian;

% BLP coefficients matrix
Bh_endBayesian=Qh_matrix'*Bh_LP+Qh_matrix_c'*Bh_VAR_Bayesian;

% Interested coefficients
Bh_endBLP_int=Bh_endBayesian(2:end,:);

% Companion matrix
Fh_varBLP = [Bh_endBLP_int';eye(n*(nL-1)) zeros(n*(nL-1),n)];

% IRFs
IRF_BLP_all(:,h+1)=reshape(bigeye*(Fh_varBLP)*bigeye'*Bzero_Bayesian,[],1);


end


end